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Abstract 

We present here a detailed description of the model of ran-driven nuclear 
transduction in living cells to be published in [9]. The mathematical model 
presented is the first to account for the active transport of molecules along the 
cytoplasmic microtubules. All parameters entering the models are thoroughly 
discussed. The simulations, carried out using the numerical scheme presented 
in [8], reproduce the behavior observed experimentally. 

1 Introduction 

All cells receive and respond to signals from their surrounding. Any external stim- 
ulus (ligand) acting on the cell plasma membrane activates several internal second 
messenger reactions that regulate virtually all aspects of cell behavior, including 
metabolism, movement, proliferation, and differentiation. As living cells are highly 
compartmentalized systems in which biochemical reactions occur in physically dis- 
tinct regions, the signal is transduced to the correct compartment to trigger the 
cellular response to the external environment. 

The cell nucleus and specifically the genomic DNA contained in it, is the target 
of many intracellular transduction pathways. Indeed the response of the cell to 
the impinging signal is obtained through the expression of specific genes. In fact, as 
protein synthesis is carried out in the ribosome, the cellular response depends as 
well on the export of RNAs out of the nucleus. 

Here, we concentrate on molecular trafficking in and out the nuclear mem- 
brane. Molecules diffuse within the cytoplasm and reach the perinuclear region. 
The translocation across the nucleocytoplasmic membrane or nuclear envelope (NE) 
may proceed through the nuclear pore complexes (NPCs) following two different 
mechanisms: passive diffusion or facilitated translocation. Only those molecules of 
mass less than 40 kDa [20] can freely diffuse through the NPCs. To permit the 
translocation of larger molecules, a system for active transport across the NPCs has 
evolved. The cargo protein equipped with a nuclear localization signal (NLS) binds 
to a nucleocytoplasmic transport receptor (NTR) karyopherin known as 'importin', 
which mediates the transport through the NE. The energy needed by facilitated 
translocation is provided by the Ran complex. The Ran protein is a small GTPase 
(for a review, see for example: [37J) cycling between two states: bound to guanosine 
triphosphate (RanGTP, active state) and to guanosine diphosphate (RanGDP, in- 
active state). The irreversible phosphorylation of RanGTP into RanGDP catalyzed 
by RanGAP in the cytoplasm maintains the Ran-cycle out of equilibrium permitting 
the accumulation of RanGTP in the nucleus [20] • This is then used to break the 
importin-cargo complex and thus permit the final release of cargo into the nucleus. 

The mechanism of facilitated translocation permits the selective and regulated 
translocation of relatively large molecules. Indeed, it is understud that the com- 
bined action of importins, Ran complex and cargoes finely regulates the action of 
transcription factors within the nucleus. 

Many mathematical models have been proposed that qualitatively reproduce the 
dynamics of intracellular signal transduction, often formulated in terms of ordinary 
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differential equations describing the time course of the molecular concentrations 
(compartmental models) [3TJ H5J [20j [3j [391 12S1 HQ] . These are based on schematic 
descriptions of signal transduction pathways which do not take into account the spa- 
tial localization of the reactions [22]. For a review on signaling networks modeling, 
see [16J. 

In these pages we discuss spatial integrated models for Ran-driven nuclear import 
of molecules incorporating diffusion and membrane transport for a large-scale model 
of living cell. For the first time, we also include a model of the active transport 
along the microtubules of the importin-cargo complex. The microtubules, together 
with the cytoplasmic filaments, constitute the cytoplasm's dynamic structure that 
maintains cell shape (cytoskeleton). The microtubules facilitates nuclear import of 
some proteins and viruses by providing a a preferential way directed towards the 
nuclear envelope (see [H] and the references therein). This is a typically spatial 
phenomena that cannot be taken into account in compartmental models. 

To the best of our knowledge, Smithed al. [45J is the only article presenting 
spatial simulations of cellular signal transduction pathways. The scarcity of spatial 
mathematical models in the literature has to be related to the lack of knowledge of 
the relevant parameters, such as local diffusion and permeability coefficients. In [45J, 
spatial and compartmental simulations are compared, giving similar answer. In our 
opinion, this has to be the case if we do not introduce any (spatial) details into the 
spatial model. In fact, if the parameters are obtained by fitting experimental data 
with compartmental simulations, spatial models may as well be less realistic. 

We shall discuss the crucial problem of parameters and pathways localization as 
we go through the various aspects of the model. 

The Ran-driven nuclear import model discussed here is presented in [9]. The 
system of non-linear equations arising from spatial modeling shall be solved using a 
new numerical technique based on Discontinuous Galerkin schemes. The details on 
the derivation and numerical properties of the method are given in [8]. 

2 The model 

The present model originates from the ODE model of Ran-driven nucleocytoplasmic 
import successively developed by Gorlich et al. (381 120] , Smith et al. [45] , Riddick 
and Macara [321 ED] , and Kopito and Elbaum [29] . 

Following [29], we keep the reaction network to the essentials. We simplify the 
set of equations and explicitly introduce the spatial component within the variables 
which shall represent the molar concentrations of the molecular species. 

Our model is a system of six semilinear parabolic PDEs set on two compartments: 
cytoplasm and nucleus. 

In the model we include, for each species, its initial concentration, molecular 
reactions, Fickian diffusion, membranes translocation conditions, and facilitated 
transport through the microtubules. 

In particular, transport through the NE and along the microtubules is only 
permitted to transport receptor complexes. 
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2.1 Reaction network and mass action law 

The reaction terms are written in terms of the Law of Mass Action (T3J. We assume 
that each biochemical reaction pathway can be decomposed into unidirectional el- 
ementary reaction and that for each elementary reaction the rate of change of the 
reactants is proportional to the product of their concentrations. The Law of Mass 
Action is in fact a mathematical model expressing the fact that the reaction de- 
pends on the number of molecular collisions and the probability that a collision has 
enough kinetic energy to initiate the reaction. The constant of proportionality is 
thus named kinetic constant. 

Experimental values of the kinetic constants are usually available from in vitro 
experimentation from purified components, and thus they do not account for envi- 
ronmental interactions, competition and localizations. In fact, one of the goals of 
spatial modeling should be to obtain the correct 'scaling relationship' of the param- 
eters by fitting to known functional effects |16j . 

Another limitation to quantitative mathematical modeling regards the rates of 
enzymatic irreversible reactions: often the literature reports the Michaelis-Menten 
(MM) kinetic parameters of the catalyzed reaction instead of the kinetic constant 
of each reaction composing it (see [331 El S3], or the review in Chapter 6 of [31]). 
Again, we have to stick to the data we are given, bearing in mind that an extra 
approximation has been introduced in the model. 

The network of reactions involved in Ran-driven nuclear import pathway are 
schematize in Figure [T] and described in Table [T] using the species names of Table [2} 



cargo 




Figure 1: Reaction scheme of Importin-mediated transport between cytoplasm and nucleus. 

We shall assume that facilitated translocation through the nucleocytoplasmic 
membrane is only permitted to molecules associated to a transport receptor (see 
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Reactions and kinetic constants for the Ran transport model 


Reaction 


loc. 


Term 


const. 


value (units) 


ref. 


Rt -? R d 


cyto 


m 1 {R t ) :- KlRg—r^ 

K *, + R t 
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20.1 s _1 

o.rjiM 


[451 


Rd -* H( 


nucl 




K d 
K t 


8.0s _1 
1.1 fiM 
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Rt + T s± T r 
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r_i(TV) := k-x T r 


fei 


0.1 (/iM s) _1 
0.3s -1 


[31 


C + T — ' T c 


cyto 


r 2 (C,T) ■- k 2 CT 




0.15 (nM s) _1 


1391 


Rt +T C -^T r + C 


nucl 


r 3 {Rt,T c ) := k 3 R t T a 


fe 3 


0.1 (/iJVT s) -1 


[ID 



Table 1: Kinetic constants for the reactions involved in the Ran-driven transport process, with 
respect to the species names in Table [2] Second order constants are measured in /itM -1 s -1 while 
first order constants are measured in s _1 . RCC1 catalyzed reaction is represented with a MM 
scheme in terms of the constant concentration cl, as in Smith e al, 2002. For a more refined scheme 
(a multistep scheme) see Gorlich et al, 2003. The MM scheme gives a good approximation as the 
intermediate steps (complex formation) are more rapid than the exchange reaction. The enzymes 
concentrations are assumed to remain constant. In particular, following[39., we set R g — 0.5/xM 
and Ci = 0.7/iM. 

below the description on membrane shuttling). 

Receptors (adapters)-mediated import necessitates RanGTP to disassociate the 
cargo from the receptor ones the complex has entered the nucleus. Thus, the trans- 
port mechanism relies on the presence of large concentrations of RanGTP in the 
nucleus in comparison to the cytoplasm (Ran gradient). The Ran complex is re- 
sponsible to maintain the RanGTP/RanGDP gradient, thus permitting cargo accu- 
mulation in the nucleus. 

Within the cell, the small GTPase Ran can bind to guanosine nucleotide phos- 
phates GDP and GTP, forming the RanGDP and RanGTP complexes, respectively. 
These reaction, described in the model by MM-kinetics, are catalyzed by two specific 
enzymes: RanGAP which is located in the cytoplasm and RanGEF (RCC1) located 
within the nucleus. Due to this cycle of reactions, a large concentration jump of 
RanGTP across the NE is readily established, with a high RanGTP concentration 
in the nucleus. 

The first cycle of cargo import begins with the formation of the cargo receptor 
complex. The receptor involved in nuclear import is made of two subunits: the 
adaptor Importin-a which binds the cargo Nuclear Localization Signal (NLS) and 
the receptor Importin-/? which permits the translocation through the NPC. Here, we 
simplify matters by considering the binding of the cargo with the whole importin 
complex (see [39J HQ] for a detailed analysis of the various import pathways). 

The second (parallel) cycle is the nuclear import of cytoplasmic RanGDP by the 
effector NTF2. Here we assume that NTF2 is always available and treat RanGDP 
as all bound to NTF2. Nuclear RanGDP, imported by the transport receptor NTF2, 
interacts with the nuclear RanGEF (RCC1) which catalyzes the exchange between 
nucleotide GDP and GTP forming the nuclear RanGTP. The complete biochemical 
scheme of this exchange is composed of a sequential set of enzymatic steps [2H1 
I2TH 135] which, following (IS] , we schematize as a single reaction. 
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The nuclear importin-cargo complex can bind the RanGTP to form two com- 
plexes, RanGTP-Imp/? and Impa-cargo. This latter complex can dissociate, create 
free Impa and cargo or interact with RanGTP to form the RanGTP-Impct complex 
and free cargo. We simplify this pathway with the single reaction 

TmC + R^C + TmR, 

written in terms of the generic transport receptor T, cargo C, and with R repre- 
senting the concentration of Ran (bound to GTP in this case). 

The free cargo, which usually is an activated transcription factor, can now bind 
DNA and activate the gene expression program. In our model, these reactions are 
neglected, and thus the cargo is allowed to accumulate within the nucleus, while the 
receptor-RanGTP complex exits the nucleus and, eventually, dissociates. 

2.2 Diffusion 

The molecular species diffusion is expressed in terms of Fick's law. All the molecular 
species diffuse with a specific diffusion coefficient. This is obtained from the formula 

D- KT 
6nr]R s ' 

in terms of the Boltzmann constant K, the temperature T, the Stokes radius R s , 
and the viscosity of the medium t] (see e.g. |5S]). 

The viscosity of the cytoplasm has been measured in the nineties for a wide range 
of molecular weights and cellular environments [44J. It was initially thought that 
the viscosity had to depend on the molecular size to account for the sieving effect 
of the cytoskeletal network and other macromolecular structures. For this reason, 
solutes diffusion was described by the translational diffusion coefficient [25] , depend- 
ing on the viscosity of the fluid-phase cytoplasm and on collisions with intracellular 
components. The former, representing the viscosity sensed by a small solute in the 
absence of interactions with macromolecules and organelles, is 1.2-1.4 times that of 
water [T8] . 

The net viscosity of the cytoplasm of Swiss 3T3 fibroblasts with respect to small 
solutes such as metabolites and nucleic acids was measured by Kao et al. [25] using 
fluorescence recovery after spot-photobleaching (FRAP) [2]. They found that mobil- 
ity was around four to fivefold slowed than in water. Contrary to expectations, later 
studies extended the validity of this result to the mobility of macro-molecules of up 
to few hundreds of kDs of mass [44j . Moreover, the same applies to the nucleus [44] 
and the mitochondrial matrix[3S], thus viscosity fourfold higher than water can be 
assumed throughout. 

A correction may be employed in the proximity of the plasma membrane as the 
translational diffusion of the fluorescent probe BCECF already used by [25] is found 
to be twofold lower near the cell membrane due to high density of proteins [47J. 

The diffusion coefficients used in the model are reported in Table [2] 
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Diffusion coefficients for Ran model 


reactant 


variable 


cliff, coef. 


value (/im 2 /s) 


reference 


RanGTP 


Rt 


d r 


22 


m 


RanGDP (with NTF2) 


Rd 


d r 


20 


|2QIS5] 


cargo 


C 


d c 


18.2 


calculated 


receptor 


T 


dt 


14 


m 


RanGTP«receptor 


T r 


dtr 


14 


m 


cargo»receptor 


T c 


d tc 


12.4 


calculated 



Table 2: Diffusion coefficients for the molecular species involved in the Ran-driven transport 
system. Values are the same for the cytoplasm and the nucleus. These data have been calculated 
assuming a cytoplasm viscosity 5 times higher than water [U] at a temperature of 20° Celsius. 
Cargo weight taken as the ERK mass, 42 kDa (for ERK1) or 44 kDa (ERK2). Data from the 
literature are not always consistent: for instance, Smith et al [45 estimates the diffusivity of Ran 
in 30^m 2 /s. 



2.3 Facilitated translocation through the nuclear envelope 

In this report, we concentrate on nuclear import by assuming that the cargo is al- 
ready in the cell periphery and that there is no exchange of substances between the 
cell and the surrounding environment through the plasma membrane. Mathemati- 
cally, no-flux conditions are specified for all species at the plasma membrane. Thus, 
given the generic species concentration u we impose: 

= o, (1) 

where n is the unit normal vector pointing outside the cell. In fact, the mecha- 
nism of signal transduction through the plasma membrane is often quite different 
from that of NE trapassing. The external signals are transduced from external mem- 
brane receptors into internal second messenger activation on the inside of the plasma 
membrane pQ. Thus, no passage of matter through the membrane is involved. 

All shuttling across the nuclear envelope takes place through the nuclear pore 
complexes (NPCs), multi-component protein structures spanning the nuclear enve- 
lope. Two mechanisms of translocation are permitted: passive diffusion or facilitated 
translocation. Due to the limited diameter of the pore lumen (about 10 nm), ions 
and small metabolites can freely diffuse through the NPCs only if their mass is less 
than ~ 40 kDa. Molecules of size greater than that of the NPCs can still shut- 
tle across the NE by facilitated translocation. The actual mechanism of facilitated 
translocation is still under research (for a review of various models, see [17]). The 
molecule must posses a specific aminoacid domain which binds to specific transport 
receptors (importins and exportins) that operate a conformational change in the 
NPC whose 'functional size' is in the order of 25 nm [32j |3]. 

Thus, we assume that facilitated translocation is allowed only to receptor com- 
plexes [2H]. 

The study presented in [2H] proves that the accumulation of cargo in the nucleus 
is not due, as previously thought, to the ability of the 'importin' receptor to cross 
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the NE is unidirectional (from the cytoplasm to the nucleus). The translocation is 
actually bidirectional, with the accumulation of cargo in the nucleus being rather 
due to the asymmetric concentration of RanGTP. 

In the absence of a detailed and realistic model of the reactions happening within 
the pore (interactions of proteins with nucleoporins, meshwork of filaments within 
the lumen of the pore, etc), a more convenient approach is to use a 'coarse-grain' 
formulation, in terms of permeabilities times concentration differences between the 
two compartments (see [13]). The flux across the nucleo-cytoplasmic boundary is 
modeled as the product of a proportionality constant (the permeability P) times the 
concentration difference across the nucleocytoplasmic boundary, as a direct conse- 
quence of the fact that the flux does not require additional energy input [35]. Thus, 
we fix the following Kedem-Katchalsky conditions at the nuclear envelope: 

d n u ^=Pu(u c -u n ), (2) 

where n is the unit vector pointing from the nucleus to the cytoplasm and the factor 
u n — u c is the concentration difference of species u across the nucleocytoplasmic 
boundary. The nuclear membrane permeability is assumed constant as experimental 
data indicate binding to the pore complex is far from saturation [3"5j H3J. The 
permeability values can be calculated from experimental values on the capacityarea, 
and number of the NPCs present on the nuclear envelope. Ribbeck and Gorlich [38J 
estimated > 100 translocation events/NPC/second and a density of NPCs of 5.1±0.2 
NPCs//im 2 in the NE of HeLa cells. 

Estimates of translocation rates of most molecular species are available from the 
literature and have already been used in a number of studies [3TJ S3J EHl ESI ED] • 
The values used in the simulations are taken from [33] and shown in Table [3j These 
were fitted by comparison of experimental data with compartmental simulations: 
more investigation on the correct values for spatial modeling is needed, for instance 
by following the approach of the homogenized 'effective' permeabilities proposed 
in HE]). 



NPC permeabilities for Ran model 


react ant 


perm, const. 


value (//m/s) 


reference (s) 


R d (with NTF2) 


Pd 


3.73 


m 


T 


Pt 


1.87 


145 




Ptr 


1.87 


145 


T c 


Ptc 


1.87 


145 



Table 3: Permeability values for molecular species involved in membrane translocation. The 
permeability is the kinetic effect of the molecular species with the nucleoporin complex (NPC) . 

The set of boundary conditions is closed by imposing continuity of the flux: 

on an 
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For molecular species which cannot pass through the nuclear membrane we just 
impose homogeneous Neumann boundary conditions on the relevant compartment. 
For instance, if the species u is confined in the cytoplasm, then 

= < 4 > 

2.4 Active transport along the microtubules 

Microtubules are known to be used by macro molecules like adenoviruses [2] and 
by intracellular organelles like endocytic vesicles [2] as preferential ways of motion. 
The adenoviruses are 90nm in diameter and thus their diffusion speed in the cyto- 
plasm is very limited: for this reason they must resort in active transport along the 
microtubules in order to reach the nucleus. 

It is by now well established that also some smaller molecules utilize motor- 
assisted transport along microtubules (MTs) (2TJ, [T|, S2] . Examples goes from mRNA 
localizing after nuclear export, to the tumor suppressor protein p53 [19] and Rb jH]. 

Many similarly sized molecules resort on passive diffusion only to reach the nu- 
cleus NE. Thus, in this respect, active transport is not essential to the transduction 
precess. Rather it must be seen as a way to improve the transduction efficiency Our 
aim is to introduce a complete model of the transduction mechanism, accounting 
for both diffusion and active transport along the micrutubules, that can be used to 
better understand the characteristics and importance of the latter mechanism. 

The microtubules, which, together with the cytoplasmic filaments, constitute 
the cytoskeleton, are usually organized into a single array centered near the nucleus. 
They are characterized by an orientation, with the minus ends located near the NE. 
Active transport along the MTs is permitted by binding to a motor protein [351 
HSj. We distinguish among two families of motor proteins: dynein, which permits 
transport from the plus end to the minus end, and kinesin, which transports in the 
opposite direction. For instance, kinesin anterograde transport along microtubules is 
known to be used by particles for covering large distances along nerve axons (see [16] 
and the references therein). 

The receptor Importin-a is able to bind to the MTs via motor proteins, thus 
permitting the active transport of NLS cargoes towards the nucleus. It is this 
phenomenon that we wish to include in our model. The speed of transport of a 
motor protein attached to an MT is of about 0.5 to 1/ims" 1 ES], but pauses 
and changes of direction of motion are often observed, suggesting that molecules 
may proceed by detaching and changing type of motor [16] . Tentative mathematical 
models to fully describe MTs effective directional transport have been proposed, for 
instance, in the References [1B1 1351 ITS] H2| 123] . 

Here, we assume that the given cargo-receptor complex can only bind to dynein, 
and thus MT transport is directed from the cell periphery to the nuclear envelope, 
and verify if the model reproduces the experimental results of Roth et al. [41J. 
We simulate active transport of Nuclear Localization Signals (NLS) along MTs by 
introducing an advection term in the species flux, which becomes: 

Ju = duVu — bu. 
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The advection field b is taken of constant modulus equal to 1 pointing towards the 
nucleus (unidirectional transport from the plus end to the minus end motored by 
dynein) . 



2.5 Partial differential equations model 

Let Q represent the cell's domain. Further, we denote by dQ its boundary (the 
plasma membrane), and by T tr the interface between cytoplasm and nucleus (the 
nuclear envelope). 

On both membranes, we fix the normal unit vector n as pointing inside the 
cytoplasm. Accordingly, the jump on T tr of a given species concentration u will be 
denoted by 

M :=u c \ rtv -u n \ rti . 

By collecting all the model's term described above, we obtain the following sys- 
tem of semilinear parabolic equations. In the cytoplasm, the species concentrations 
obey: 

d r AR t - m x (Rt) - n{Rt, T) + r_i(T r ), 



dt 
dRd 

dt 
8T r 

~dt 
dC 
~dt 
dT 
~dt 
dTc 
dt 



d r AR d + m^Rt), 

dtrATr + niRt^-r^Tr), 

d c AC-r 2 (C,T), 

d t AT - r x {R t , T) + r_x(T r ) - r 2 (C, T), 
d tc AT c -V(bT c )+r 2 (C,T), 



Notice the advective term modeling transport of the receptor-cargo complex which 
appears in the last equation above. 

As we assume that no matter is entering or exiting the cell through the plasma 
membrane, the homogeneous Neumann boundary condition Q is imposed to all 
species. 
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In the nucleus, we have the following system of reaction-diffusion equations: 
( dR t 



dt 
dR d 

dt 
dT r 

~dt 
dC 

~dt 
dT 

~dt 
dT, 
dt 



d r AR t + m 2 (R d ) - n(R t ,T) + r^(T r ) - r 3 (R t} T c 
d r AR d - m 2 {Rd), 

d tr AT r + n(R t , T) - r_i(T r ) + r 3 (R t , T c ), 

d c AC + r 3 (Rt,T e ), 

d t AT — ri(R t , T) + r_i(T r ), 

d tc AT c + V(hT c )-r 3 (R t ,T c ), 



The two systems of equations above are coupled through the appropriate trans- 
mission conditions on r tr : 



( A dR T\ 

«r 9n |r tr - 


0, 




A 9R d\ 

a r-Q^\r tI - 


A 9R d\ 




d aT -L - 




Ptr{T r 


j dC c - n 1 

«c 9n |r tr - 






"*a^|r tf - 


"fair Ir* - 




#7 9Tcc lr, - 

I u * c an irtr — 


»tc lr tr - 


Ptc{T c 



These equations needs to be provided with initial conditions which are specified 
below in Section HI 



3 Numerical approximation of the mathematical 
model 

We name the two subdomains of Q representing the cytoplasm and the nucleus as 
Q\ and Q 2 , respectively. Thus, Q = Qi U 5722 U r tr , where T tr := £l\ fl fl 2 - 

In order to ease notation in the presentation of the numerical method, we rewrite 
the mathematical model described in the previous section in vector form. 

3.1 The pde problem in vector form 

Let us denote the number of unknown concentrations by n, and collect all concen- 
trations variables in the vector function 

u :=(«!,... , u n ) T i^Ufia^ M n . 
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We shall look for solutions belonging to the space 7i l := [H l (Vt\ U fi 2 )] n a t any 
time t G [0, T], with T representing some final time. 

We introduce the diagonal matrix U = diag(«i, . . . ,u n ) and the gradient Vu : 
fiiUfia^ M nxd given by Vu := (Vui, . . . , V« n ) T , with Vuj(x) G R d , i = 1, . . . , n, 
x G Qi U f2 2 . Further, for a tensor Q : Hi U ^ -> IR nxd , with rows Qi, i — 1, . . . ,n, 
we define its divergence V • Q : fii U fi 2 ->• K n by V • Q := (V • Q u . . . , V • <2„) T . 

We collect all reaction terms in the vector field f(u). Similarly, the advection 
and diffusion coefficients are collected in the tensors B G [C 1 (0, T; 0\r tr )] Tlx ' i , whose 
rows are denoted by B^, i = 1, . . . , n, and A G [C(0, T; Qi U f2 2 )] nxn diagonal, with 
A = diag(di, d 2 , . . . , d n ), where di > 0, i = 1, . . . , n are the diffusivity of the various 
species. Notice that Bi will be the zero vector if the z-th species is not transported. 
In particular, we assume that every Bi is continuous and supported on an open 
(eventually empty) subset of Q\ (the MTs are located in the cytoplasm). The 
implication of this assumption which is used below in the definition of the method 
is that B = on the boundaries of both Q\ and f2 2 . Finally, let uo G [L 2 (fi)] n 
represent the initial conditions. 

We look for solutions u G L 2 (0,T;7i 1 ) of the following initial and boundary 
value problem: 

' u t = V • (AVu - UB) + f (u) in [0, T] x U fi 2 ), 

u(0, x) = u (x) on {0} x Q, 

< (AVu)n = on dQ, (5) 

(AVu)n| ni = P (u|n 2 - u|nj on T tr , 
k (v4Vu)n|n 2 = P (u| ni - u|n 2 ) on r tr , 

where n|n. denotes the unit outward normal vector to r tr from Qj, j = 1,2, and 
P := diag(pi, . . . ,p n ) the diagonal matrix of permeabilities p^. 

More general boundary conditions, including mixed-type conditions, and trans- 
mission conditions, including nonlinear conditions, are discussed in [8]. 

3.2 The interior penalty discontinuous Galerkin method 

In recent years, there has been an increasing interest in discontinuous Galerkin 
finite element methods (DGFEM)[12J. The renewed interest in these methods is 
due to their very good stability properties when used to approximate solutions to 
convection-dominated convection-diffusion problems [11] , as well as due to the great 
flexibility in grid design they offer. Moreover, they naturally embed good local 
conservation properties of the state variable, which can be advantageous for time- 
dependent problems. 

Regarding the application of DGFEMs to the present problem, we add that the 
weak imposition of boundary conditions typical of discontinuous methods permits a 
very natural imposition of the transmission conditions. The good stability properties 
of the method on convection-dominated problems are relevant to the case of cargos 
representing large viruses, not considered in the simulations below. Indeed, the 
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diffusion of macromolecular structures is impeded by the cytoskeleton (the effective 
viscosity is up to a hundred times higher than that of water), and active transport 
along the MTs becomes the dominating mean of migration towards the nucleus [7] . 

The method used here, which is detailed in [8], is a generalization of the sym- 
metric version of the interior penalty discontinuous Galerkin (IPDG) method. This 
was introduced in [23] for the solution of second-order partial differential equations 
with nonnegative characteristic form and later extended to the solution of semilinear 
time-dependent problems in |3Uj . 

The finite element method is introduced by defining a shape-regular subdivision 
T of fl± U f^2 into disjoint open elements k G T. We assume that the elements of T 
'belong' to either subdomain, so that subdivisions of f2i and Q2 are automatically 
introduced. Abusing notation, we use the simbol r tr to denote the union of elemental 
faces belonging to both subdivisions. Further we decompose the subdivision skeleton 
T := U Ke r<9K into three disjoint subsets 

r = dvt u r int u r tr , 

where T int := T\{d£l U T tr ). 

We assume that the subdivision T is constructed via mappings F K , where F K : 
k — > k are smooth maps with non-singular Jacobian, and k is the reference d- 
dimensional simplex or the reference <i-dimensional (hyper)cube; the maps are as- 
sumed to be constructed so as to ensure that the union of the closures of the elements 
kG T forms a covering of the closure of Q, i.e., Cl = U Ke rk. 

For a nonnegative integer m, we denote by V m {k), the set of polynomials of 
total degree at most m if k is the reference simplex, or the set of all tensor-product 
polynomials on k of degree k in each variable, if k is the reference hypercube. We 
consider the /ip-discontinuous finite element space 

S := {v G L 2 (n) : v\ K o F K E V mK (k), « e T}. (6) 

Next, we introduce some trace operators. Let k + , k~ be two (generic) elements 
sharing an edge e := 8k + H C r int U r tr . Define the outward normal unit 
vectors n + and n~ on e corresponding to dn + and 8k~, respectively. For functions 
q : Q —>■ W 1 and Q : Q —>■ W axd that may be discontinuous across T, we define 
the following quantities. For q+ := q\ eC d K +, q~ := qUcd*-, and Q + := Q| eC a K +, 
:= Q\ecd K -, we set 

{q} :=^(q + + q-), {Q} := ^(Q + + Q"), 

and 

Iq| := q + ® n+ + q~ ® n", [Q] := Q+n+ + Q n", 

where <S> denotes the standard tensor product operator, whereby q ® w = qw T . If 
e G Ok Pi Tg, these definitions are modified as follows 

{q} := q + , {Q} := Q + , M := Q + ® n, [Q] := Q + n. 



14 



Finally, we introduce the mesh quantities h : Q — > K, m : Q — > R, by h(x) = 
dianiK, m(x) = m K , if i e k, and h(x) = {h}, if x G T, m(x) = {m}, if x G T, 
respectively. 

In order to define the IPDG finite element method, we further introduce £ := 
diag(<7i, . . . , a n ) the diagonal matrix containing the discontinuity-penalisation pa- 
rameters &i : r\r tr , and denote B := 1/2 diag(|£>i • n|, . . . , \B n • n|). 

The IPDG-in-space method for the system (|5| reads: find u/j G L 2 (0,T; \S] n ) 
such that 

((u ft ) t ,v fc )+B(u h ,v h ) = (f(u fc ),v fc >, for all v fe G L 2 (0, T; [5]"), (7) 
where (•, •) is the L 2 scalar product, the bilinear form and B(-,-) is defined by 

B(u h ,v h ) ! (AVu h - U h B) : Vv ft + f P[u h ] : [v h ] 

- / ({iVu, - LThB} : [v A ] + { AVv,} : |uj - (E + B) fuj : [vj) . 

•J Tint 

(8) 

A fully discrete formulation is obtained from ^ by using any standard time step- 
ping method like, for instance any Runge-Kutta time stepping. Full details on the 
derivation and analysis of the method can be found in [8]. 



4 Numerical simulations 

We use the numerical method described above to simulate our model of cargo import. 
The goal is to verify the experimental results presented by Roth et al. [H]. Two 
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Table 4: Initial concentrations of constituents of the Ran system used in the model by Riddick et 
al. [35]. All other initial concentrations are set to zero. The concentrations of the enzymes RCC1 
and RanGAP are assumed to be constant, as described in Section [27T] 

experiments are presented in [UJ: in the first, the NLS cargo known to bind to the 
MTs is activated in a cell with intact cytoskeleton, while in the second experiment 
the cell is treated with the MT-depolymerizing agent nocodazole (NCZ) to create 
a MT-less environment. In silico we can easily turn on and off the advection term 
and compare cargo accumulation in the nucleus. Notice that this corresponds to 
permitting or not the association to the MTs, while depolymerization cancels the 
MTs thus changing as well the cell environment. The full consequences on the 
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cell metabolism under depolymerization are unknown, and thus comparison with in 
silico experimentation is crucial. 

The initial conditions are the concentrations of the single molecular species when 
cell is at rest (i.e. in the absence of external stimuli). Initial concentrations used in 
the simulation are reported in Table |4| In particular we assume that the activated 
cargo is initially concentrated in a peripheral zone of the cytoplasm, as shown in 
the left plot in Figure [2] In such a situation, the concentrations of the complexes 
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Fi gure 2: NLS cargo concentration in the cell at the initial time (left) and after 17 seconds (right). 

are zero (i.e. no complex is formed before the stimulus activates the cargo). The 
experimentally observed RanGTP accumulation of RanGTP in the nucleus forms in 
the initial phase of simulation. 

The accumulation in time of cargo mass in the nucleus is shown in Figure [3j The 
change in accumulation rate is evident, and confirms the behavior experimentally 
obtained in |41j. 
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Figure 3: Accumulation of the cargo in the nucleus. 

A snapshot showing the simulated cargo concentration in a spherical cell after 17 
seconds is shown in Figure [2] (right), where the accumulation jump across the NE is 
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clearly visible. Figure [4] shows the concentration after the same time lapse, this time 
obtained with a two dimensional computation to better appreciate the variation of 
concentration along the cell. Notice in particular that the concentration inside the 
nucleus is, in presence MTs active transport, more than twofold. 




Fi gure 4: NLS cargo concentration in the cell after 17 second. With advection along the MTs(left) 
and without (right). 
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